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EPIDEMIC SIZE IN THE SIS MODEL OF ENDEMIC INFEC- 
TIONS DAVID A. KESSLER,* Bar-Ilan University 
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Abstract 

We study the Susceptible-Infected-Susceptible model of the spread of an 
endemic infection. We calculate an exact expression for the mean number of 
transmissions for all values of the population and the infectivity. We derive the 
large- A'^ asymptotic behavior for the infectivitiy below, above, and in the critical 
region. We obtain an analytical expression for the probability distribution of 
the number of transmissions, n, in the critical region. We show that this 
distribution has a n~^^'^ singularity for small n and decays exponentially for 
large n. The exponent decreases with the distance from threshold, diverging 
to infinity far below and approaching zero far above. 
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1. Overview 

The Susceptible-Infected-Susceptible (SIS) model of Weiss and Dishon [15 is one 
of the simplest model of endemic infections. The model describes the evolution of an 
infection in a fixed population, with no restriction on the possibility of reinfection of 
a previously infected and now recovered individual. This is contrast to the venerable 
Susceptible-Infection-Recovered (SIR) model [8^, where reinfection is not permitted. 
Both models exhibit a threshold value of the infectivity, below which the infection im- 
mediately dies out. Below threshold, then, where only a tiny fraction of the population 
is impacted, the two models have essentially equivalent statistical properties. Above 
threshold, in the SIR model the infection is self-limiting, since in a fixed population 
the number of potential new victims, the susceptible pool, is monotonically decreasing 
in size. The SIS model, on the other hand, describes an endemic infection which can 
(above threshold) persist indefinitely, at least at the deterministic level. Thus, the 
statistics of infection size in the two models above threshold are very different. 

The statistics of the mean time to extinction in the SIS model have been much 
studied, starting with the original paper of Weiss and Dishon |15| and most recently 
by Doering, Sargsyan and Sander ^ . The latter paper investigates the large population 
limit of the mean extinction time. This goes from a logarithmic dependence on the 
population size, N, below threshold, to a V-ZV dependence exactly at threshold to 
an exponential dependence above. In this paper we will focus on the mean number 
of transmissions till extinction. This is a more pertinent method of characterizing 
the epidemic and the threshold transition. We shall derive an exact formula for this 
quantity for general N and infectivity, and then examine its large N asymptotics. As 
we shall see, above threshold the mean epidemic size is directly related to the mean 
epidemic duration. At and below threshold, though, these quantities are quite different. 
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Furthermore, the number of infection events is directly relevant when considering the 
probability of a mutation of the pathogen, as mutations are most probable during the 
exponential growth phase following a new infection [5]. These mutations are implicated 
in the conversion of a sub-threshold weakly transmittable pathogen into a super- 
threshold variety capable of inducing a major epidemic. 

Of particular interest will be the critical regime separating the above and below 
threshold cases. As already noted by Nasell [13], for a range of infectivities of width 
1/VN around threshold, there is a crossover region that interpolates between the above 
and below threshold cases. The existence of a large N scaling theory in this region 
was recently proven by Dolgoarshinnykh and Lalley E]. We shall see this crossover 
region and its characteristic scaling arising naturally from our general result for the 
mean infection size. 

After this treatment, dealing exclusively with the perhaps most biologically relevant 
case of a single initial infection, we extend our results to an arbitrary number of 
initial infections, again deriving an exact formula and then examining the large N 
asymptotics. In the crossover regime, we will have to distinguish the cases when the 
number of initial infections in small, comparable to, or much larger than ^/N- 

From looking only at the mean number of infections, we move on to consider the 
entire probability distribution for the number of infections. We first briefly discuss the 
above and below threshold cases, and then focus in on the critical threshold regime. In 
the particular case of exactly at threshold, the entire probability distribution for the 
appropriate scaling variable (the number of infections divided by N) can be explicitly 
displayed. In general, we can express the probability distribution as an inverse Laplace 
transform. This is sufficient to calculate the limiting behavior of the distribution for 
small and large epidemics, and to recover our expression for the mean in the critical 
regime. We then conclude with a few observations. 



2. Preliminaries 

We begin with a description of the SIS model. The individuals in the population 
are divided into two subclasses: the susceptible pool, of size S, and the infected (and 
infectious) class, of size /, with A^ = S + I. The disease is transmitted from an infected 
individual to a susceptible one with rate a/N, so that 

(5,/) "^-^/^ (5-1,7+1). 
Infected individuals recover with a rate /3, reverting back to susceptibles: 

(5,7)^(5+1,1-1). 

Of primary interest is the case where initially 5 = A^— 1,7=1, so that the outbreak is 
sparked by a single infected individual. The outbreak terminates when the last infected 
individual recovers, and 7 returns to 0. 

This stochastic process is traditionally approximated (for large populations) by the 
rate equations 

5 = -^57 + /37 
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Using the conservation of A^, we get 



I ^ {a- 13)1 



a 



N 



which is a logistic-type equation. We see that there is a transitition at Rq = a/f3 = 1, 
where i?o is equal to the mean number of primary infections caused in a large population 
of susceptibles by an infected individual. It is clear that if i?o < 1, the 7=0 state 
is stable, whereas for Rq > I the rate equation predicts a stable equilibrium state 
at 5* = {f}/a)N, — (I — f3/a)N. Thus at the classical level, Rq — 1 marks the 
threshold between an infection that becomes endemic and those that fail to spread. 



Already in the original Wciss-Dishon paper \i5\, an exact expression for the mean 
time to extinction, starting from the completely infected state, was derived. The 
generalization of this to an arbitrary number of initial infected individual was given 
in Leigh 111] and rediscovered by Doering, et al. [1]. However, the mean number of 
infections is the quantity of primary interest in an infection model. We can focus in on 
this quantity if we eliminate time, considering only the transitions between states. We 
characterize the system by the number, T, of transitions the system has undergone. 
In each transition the number of infected individuals either rises or falls by one, so 
that / undergoes a kind of random walk. The probability of an upward transition is 
p+ = RoS/{RaS + N) = Ro{N - I)/{Rq{N - I) + N), whereas the probabiHty of a 
downward transition is = 1 — p^. These probabilities are unequal and depend on 
J, so that the walk is biased, with a "space" -dependent drift. (From here on, we will 
refer to T as Time, with the lower case word "time" retaining its usual meaning, and 
trust this will not lead to confusion). It is easy to see that at the point of extinction, 
the total number of infections, including the initial Uo infected individuals, is just 
n = {Text + no) /I. The number of induced infections is of course Uo smaller. Since the 
results of Weiss and Dishon, Leigh and Sander, et al. for the mean time to extinction 
apply to a general one-step random walk, we can apply them directly to calculate the 
mean number of infections. Specializing to the case where we initially have exactly one 
infected, we have for the mean extinction Time, ri 



Here we have indicated explicitly the dependence of the transition probabilities on / 
via a subscript. Plugging in these probabilities, we find 



3. Mean Number of Infections 





Reordering the sum, we can rewrite this as 



Tl = 




Ron\ 1 ( N 



N ) n\ \R^ 
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We can do better, since the second term in the sum in the same as the first, except for 
the last index, so 



1 / N 



N-l 



SO that 



We recognize the sum as the first N terms of the Taylor expansion of the exponential 
^N/Ro^ The behavior of the sum depends on whether Rq is above or below 1. This 
follows from the fact that the terms in the Taylor expansion of increase until n — x, 
and then decrease. For large x, in fact, the behavior of the terms with n is a Gaussian 
peaked at n = a;. The behavior of the sum is then determined by whether the last term 
of the sum at rt = — 1 occurs before or after the peak at rt = N/Rq, i.e. whether Rq 
is above or below 1. 

For i?o above 1, the summed terms extend past the peak, which dominates the 
sum, and so, up to exponentially small corrections, the sum is just the exponential. 
Furthermore the prefactor can be approximated via Stirling's formula, giving 

_^ _ V^^^N{lniRo) + l/Ro-l) (X) 

Ro 

Thus, as expected the mean number of infected cases grows exponentially large with 
TV, with the exponent going to as Rq approaches 1. Furthermore, the exponent is the 
same as for the mean first passage time (here actual time) as calculated in Ref . |1] , and 
is equal to the action for the semiclassical extinction trajectory [6,|9j. This is because 
above threshold, the system remains an exponentially long time in the classically stable 
state. We plot n versus Rq in Fig. [T] together with the large-TV asymptotic formula, 
Eq. 0. We see that the mean number of infections quickly grows to astronomical 
proportions as Rq increases away from 1. To see the approach to the large-iV result, 
we show in the inset the ratio of the exact results for N = 100 and 400 to the large- iV 
asymptotic formula. We see that the approximate formula works excellently except in 
the vicinity of the transition point _Ro = li and improves with increasing A'^. 

For Rq below 1, the sum is cut off while the terms are still increasing with n, and so 
the largest terms in the sum are the last ones, which approximate a geometric series 

m«V4 = 



fc=0 " 



which of course is the same as in the SIR model, since the number of infected persons 
is so small that no one gets a multiple infection. This infinite N answer is compared 
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Figure 1: Supercritical Regime (Rq > 1); Mean number of infections, n versus the infectivity 
parameter _Ro, for the case of starting for one infected individual. Also shown is the large A'' 
asymptotic formula, Eq. ([T|. The cases of population A'^ = 100 is shown. Inset: Ratio of the 
exact mean number of infections, n to that given by the large N asymptotic formula, Eq. ([l]), 
versus the infectivity parameter Rq. The cases of population N = 100, 400 are shown. 

to the finite N results in Fig. |2] where we see that it works well as long as we are 
sufficiently below Rq = 1, and the range of agreement increases with N. As opposed 
to the supercritical case, here there is in general no simple relation between the mean 
number of infections and the mean (actual) time to extinction, which is given by 



1 



N 



77 



(^-1)! 

{N-jy. 



fc=i 



1 



ln(l - Ro) 



This is of course due to the fact that in the subcritical case the number of infections 
in not strongly peaked about some value as it was in the supercritical case. 

For Rq near 1, i?o = 1 + 5N~^/'^, there is a transition region. The dominant terms 
in the sum are again the largest, which have a Gaussian character, with a maximum 
at N/Rq « N. 
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Figure 2: Subcritical Regime {Ro < 1): Mean number of infections, n, versus the infectivity 
parameter Ro, when starting with one infected individual, together with the prediction for an 
infinite population. The cases of population N = 100, 400, 1600 are shown. 

This clearly reproduces the sub- and supercritical results in the limit of large negative 
and positive S, respectively. This formula is plotted in Fig. [3] along with data for 
N = 100,400,6400. We see that finite N data converge to the infinite N prediction, 
with the finite N effects larger at larger S. 

4. Mean Number of Infections, General Initial Condition 

These results are easily generalized to the case of Uq initial infected individuals, 
again starting from the corresponding mean first passage Time. We get 



ni 



k=2 



0-fc 



where 



CTfc 



~N 



N-k 



N-k 



1 



^ j! 



{N-ky. X" 

Again these results are instructive in the various limits. For the above threshold case, 
R]f~^ni so that 



1 



1 
i?0 



1 



1 



, 1 — i?n 



Ro-l 



,7V(ln(flo) + l/flo-l) 



(3) 



The prefactor is recognized as the probability of a biased random walk starting at Uo 
to survive to infinity. Thus, the mean number of infections is the mean number of 



Short title 



7 




Figure 3: Critical Regime (_Ro = 1 + 5/\/N): Scaled mean number of infections, n/\/N 
versus scaled infectivity parameter 5, when starting with one infected individual. Also shown 
is the large N asymptotic result, Eq. ([2|. The cases of population A'' — 100,400,6400 are 
shown. 

infections starting in the macroscopically infected metastable state times the proba- 
bihty of surviving long enough to reach this state. The expression in brackets, the 
mean number of infections starting in the metastable state, is itself simply 

related to the mean first passage time for this initial state calculated in Ref. |1] . If one 
accounts for the average time for a transition in the metastable state, 2(i?o — l)/-^-Ro, 
one can easily obtain from the above the average Time to extinction, which is twice 
the average number of infections. This is because the overwhelming majority of the 
time is spent in the vicinity of the metastable state. In Fig. |4] we present the exact 
results for for the case N — 100, Rq = 3 together with our approximation, Eq. 
([3]). We see the agreement is quite satisfactory. 
Below threshold. 



00 r / 7 N -II 2 



1 - (1 - k/N)Ro 



Thus, for Uo N all the cr's are all approximately equal to fii, so that « Uofii. 
This is clear, as the individual seeds act essentially independently, since they impact 
an infinitesimal fraction of the total population. This is in sharp contrast to the above 
threshold case, where h„^ converges exponentially quickly over an 0(1) range of tIq. 
For larger Hq of order N, we get 

For the subcritical case, it is also interesting to consider the number of induced infec- 
tions, since here most of the infections are just those of the initial state. For small Uq 
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Figure 4: Exact calculation of n as a function of Uo for the supercritical case, _Ro = 3 together 
with the analytical approximation, Eq. ([3]|. A'^ = 100. 

the number of induced infections is approximately noi?o/(l — ^o), again proportional 
to Uo- For Uo = N, the number of induced infections is (A^/i?o) ln(l/(l — Rq)) — N, 
which for small Rq is approximately NRo/2. Thus the interference between different 
initial seeds reduces the number of induced infections roughly by half in this case. The 
interference effect is of course even more dramatic for larger Rq. As Rq approaches 
unity, the number of induced infections diverges only logarithmically for Uo = N, as 
opposed to the 1/(1 — Rq) divergence for small Uq. Of course, for Rq even larger, 
in the supercritical regime, as we have seen, the interference effect is almost total, 
as increasing Uq beyond 10 or so has essentially no effect. The subcritical case is 
demonstrated in Fig. |5]for the case Rq — 0.3. 

In the critical regime, things are of course a bit more complicated. Again, we first 
compute the ak- 



Af-fc j-1 

j=o e=Q 



fc + i 

N 



,(5-fc/V]V)V2 




We now have to integrate this with respect to fc, whose typical scale is 0(v7V): 



fir,. = iV 



l + crf(^(5-x) 



(5) 
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Figure 5: Exact calculation of n as a function of Uo for the subcritical case, Ro ~ 0.3 
together with the analytical approximation, Eq. Q. Also shown is the average number of 
induced infections, n„„ — Uo-N = 100. 



One immediate result is that n„^ starts out as 0{\/ N) for small Uq of order unity, 
but for Uo of order 0{^/N), the average number of infections rises to 0{N). In Fig. 
[6] we present the exact results for n„^ in the critical region for N = 400 versus our 
scaling prediction, Eq. (|5|. We see that the scaling results are perfect for the exactly 
critical case, and the further we are from criticality, the larger the finite N effects are. 
Furthermore, the finite N effects are larger for positive 5 than for negative. Also, the 
larger the initial infection size, the larger the finite N effects. 

To better understand our scaling formula, we consider in turn the cases 5 large and 
negative, 5 of order unity, and 5 large and positive. In the former case, we the argument 
of the erf is large and negative and so 



dk- 



1 



= TV In 1 



\5\ + k/^N V ■ l^h 

Thus, n„^ crosses over from a linear behavior at Uo of order unity, to logarithmic 
growth when Uo of order ^/N . For (5 = of order unity, it is more useful to integrate 
first with respect to k and then do the integral over j. This gives the formula 

dje 



N 



Again, for Uq small compared to \N, the answer is proportional to Uo- For large Uo, 
this can be approximated as follows: 



hm N 



N 



(6) 
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Figure 6: Scaled average infection size n/N as a function of the scaled initial number 
of infected individuals, Uo/VN for iV = 400 and 5 = -2,-1,0,1 {Ro = 0.9,0.95,1,1.05 
respectively). 

where A{S) is given be 

tF ,2 /T , f°^' cosh fc(5 - 1 



A is a monotonically increasing function of 5, with A(0) — 0. For large negative 6, 
A{d) ~ — ln((5) — 7/2 — (ln2)/2, reproducing our previous result. For large positive S, 
A grows quickly, A{S) « V2ne^ Z^. Thus, for all S, n grows logarithmically in Uo for 
no/V^ sufficiently large. However, since A grows so rapidly with S, for large positive 
S this behavior is not readily visible in practice, since Hq can be no bigger than N. In 
Fig. |7] we show h„^ for different S's. We see that the large rio approximation works 
well for no/\/N larger than 1 or so. For large positive 5, the argument of the erf is 
large and positive, and so the behavior for fixed Uo/y/N is most relevant. Then, 



The problem with this expression is that it is not at all accurate until 6 is fairly large, 
around 6 or so. For such large ^'s, the concept of a critical region does not apply 
until really large N's. For example the next order correction to the action is of order 
S^/N^^"^, which is only small for N ^ S^. We can see this in Fig. [s] where we examine 
the convergence of the finite- results to the critical scaling result for 6 — 3. Only for 
A'^ = 160, 00 is h approaching its limiting scaling form. This form is well approximated 
by our large Uq formula, Eq. ^ for rtol^fN > 2, and by the large 5 formula, Eq. ^ 
for Uo/y/N < 1/2. 
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Figure 7: Calculation of scaling form of scaled mean epidemic size n/N as a function of 
scaled initial epidemic size Uoj^fN in the critical regime for intermediate 5 = —2, —1, 0, 1, 2. 
Also plotted is the large rio analytic approximation, Eq. ||6|. 

5. Distribution of Number of Infections 

We now turn to investigate the distribution of the total number of infections, 
returning again to the case of a single initial infected person, Uq = I. For the sub- and 
supercritical cases, these are fairly simple. In the subcritical case, almost surely the 
infection goes extinct before the percentage of infecteds, I/N, is significant. In this 
case, the up transition probabilities are essentially constant: — Ro{N — I) / {Ro{N — 
I) + N) w i?o, PjT = 1 — « 1 — -Rq- Therefore, the random walk reduces to that of 
a constant leftward bias, with P„ falling exponentially with n: 



P{n) = 



T>n- 



(1 + i?o)2"-l 

1 (4flo)"~^ 
V^(l + i?o)2"- 



2n-2 
n- 1 



2n - 2 
n 



{n > 1) 



(8) 



Above threshold, the system spends an exponentially long time in the metastable 
state. Thus Pin) for macroscopic n's is an exponential distribution, the waiting time 
distribution for the decay of the metastable state. For small epidemics, where the 
space-dependent drift is not yet relevant, the system can again be approximated by a 
random walk with constant bias, this time to the right. Thus for n <^ N, P{n) is given 
by Eq. For larger N, P{n) crosses over to a pure exponential decay, normalized to 
1 — I/Rq, the probability of the infection surviving to macroscopic size. This behavior 
is exhibited in Fig. |9] It is important to note the difference between this behavior and 
that exhibited above threshold in the SIR model [TH|T5J|TD]. There the distribution has 
a second peak (in addition to the one at the origin) at the number of infections predicted 
by the deterministic dynamics. In the SIS model, the total number of infections above 
threshold predicted by the deterministic dynamics is infinite. Rather the behavior for 
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Figure 8: Scaling average epidemic size n/N as a function of scaled initial epidemic size Uo 
for N = 10", 1.6 • lO'^ for the above critical case (5 = 3. Also shown in the analytic critical 
scaling form, Eq. ([5|, and its large rio (Eq. ([6|) and large 6 (Eq. ([7|) limits. 



major epidemics is the pure exponential waiting-time distribution, with its peak at the 
origin. 

We now turn to investigate the behavior in the crossover regime, Rq = 1 + . 
The essential simplification here is in the transition probabilities, which in this regime 
can be approximated by 

1 / 5 

Pi = - =F ± — ^ 

^ 2 4iV AyfN 



The critical regime is characterized by the scaling / ~ 0{N), so that the bias is small, 
of order . For relatively small Times, (n ^ \/N) the bias is irrelevant, and the 

problem reduces to the unbiased random walk, given by substituting i?o = 1 in Eq. 
^ above: 



P{n) 



1 



22n-l 
1 



'/2n 


:0 






\ n 









(n> 1) 



(9) 



We now study how for larger Times the bias, resulting from the reduction of the 
susceptible pool with increasing / and the small deviation from criticality, modifies 
this answer. 

As the bias is very weak, however, and only effective at large Times, we are jus- 
tified in passing to a Fokker-Planck description for the the probability distribution 
K{I, T; lo), for /, given that there were lo infected individuals at T = 0. For a typical 
major epidemic, T at extinction is of order N, this despite the fact that n is of order 
-\/iV, since the probability of a major epidemic is of order 1/^/N. We thus define 
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Figure 9: Probability distribution P{n) for total size of epidemic in the supercritical case, 
Ro = 1.2, A'' = 300. Left: Behavior for n < N , together with the constant bias approximation, 
Eq. ([8|. The beginning of the crossover to exponential behavior is visible toward the end. 
Right: Behavior for n on the scale of the average epidemic size, together with an exponential 
fit. The normalization is seen to be 5.510~^ • 3.710~^ = 0.2, consistent with our expectation 
of 1 - 1/1.2 = 0.17. 



t = T/{2N) and x = I/VN and consider K{x, t; Xo): 
d 



dt 



92 d dK 



(10) 



with the initial condition K(x, t; Xq) — d{x — Xo) and the absorbing boundary condition 
if (0, t; Xo) — 0. In terms of K, the probability distribution for the epidemic size, P{n), 
is given by 

P{n) = K{x,n/N\Xo) 



dxdxn 



(11) 



since Xo — Iq/VN — l/^/N. This critical regime distribution function, obtained from 
a numerical solution of Eqs. ( 10 1 and (111, is presented in Fig. 10 for 6 = —1, 0, and 1. 
We see that for small n all three curves collapse into a universal power law. For large 
n the distributions fall off rapidly, with the speed of falloff decreasing with increasing 
5. 

To make more progress, we recognize Eq. (10 1 as the imaginary-time Schrodinger 
Eq. for a harmonic oscillator potential, modulo a Gaussian similarity transformation. 
Defining 



K{x, t; Xo) = e^—o)s/2-i.^-.t)/i+t/2Q^^^ ^. 



we have 



G=G" -\{x-SfG 
4 



with G{x, 0; Xq) — S(x — Xq), G(0, t; xq) = in terms of which 

92 



P{n) = iV-3/2en/2W 



dxdxo 



G{x, n/N; Xo) 



The only complication is the presence of the absorbing boundary condition at a; = 0, 
which breaks the reflection symmetry of the potential around x ^ S. Exactly at 
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Figure 10: N^^'^ P{n), the scaled probability density of epidemics of size n in the critical 
regime for the cases 5 = — 1, 0, and 1, as a function of scaled epidemic size n/N . 



theshold, 5 = 0, however, this is not a problem, as the boundary condition can be 
enforced by the method of images. The solution is 

G{x, Xo, t) = A{t) sinh {2c{t)xx„Xt)) e-'=(*)(=^'-^™(*)') 



where 



•^m (^) 

A{t) 



1 



cotht 



cosh t 
1 



sinh^i/^te-Ttanht 



This gives the probability distribution 
P{n) = , ^ 



e"/2^sinh-3/2(n/iV) 



(12) 



This clearly reproduces the expected behavior, Eq. Q, for 1 <C n ^ A^, and then 
decays exponentially for n ~ N . This result is shown in Fig. [TT] together with data 
for N = 25 and 100. Even for these small TV's, the agreement is excellent, exact at the 
smallest n's, where the discreteness of n factors in. 
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Knowing P(n) gives us another way to calculate the mean epidemic size, n: 

nP{n) = \ f dx- 



2N , 

dx- 



Jo Vl 



e 



sm i(e -)|^ 



ttN 



which of course agrees with our previous result, Eq. ([2]), specialized to S — 0. 

We now return to the distribution of n for general 6. Formally, we can decompose 
G{x,t; Xo) into a sum over eigenfunctions of the Shroedinger operator: 

as follows: 



G{X, Xo, t) = ^ (f>n{x)(j>n{Xo)e~ 



■E„t 



If we define the Green's function G(x, x' , E) as usual by 

(j>n{x)cj)n{x') 



G{x,x',E)^Y. 



En — E 
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then G{x^ x\ —E) is the Laplace transform of G{x, x\ t) with respect to time: 

POO 

G{x,x',E)^ dte^*G{x,x',t) 
Jo 

and G{x,x',E) satisfies 

HG-EG^ 5{x - x') 
We can recover G{t) from G{E) by an inverse Laplace transform 

/'7+ioo 



1 PJ + IOO 

G{x, x', t) = — dEe-^'G{x, x' , E) 



where 7 lies to the left of all the poles of G. Defining P{E) as: 

92 



P{E) 



dxdx' 



-G{x,x',E) 



x—x'—O 

we get an expression for our desired probability distribution P{n)\ 



P{n) = 



1 



27ri7V3/2 



7 + 



7 — 200 



,-{E-l/2)n/N 



P{E) 



Thus all we need to do is calculated the Green's function G{x, x' , E). There is a 
nice formula relating G{x,x',E) to G^{x,x',E), the Green's function of the problem 
without the wall. The Green's function for the system with a wall on the left-hand 
side of the system at a; = a is [7] 



G{x,x',E) = G°{x,x',E) - 



G"{x,a,E)G°{x\a,E) 
GO{a,a,E) 



Denoting f{x,E) as the solution of Hf ~ Ef which decays as x ^ +00, and g the 
solution which decays as a; — s- — cx), then in general 

rO(r r' E) - ■^(^>)g(^<) 



where Wr[f,g] is the Wronskian. Then, using the fact that the FFr[/, g] is constant, 
P{E) 



x—x'—O 



1 



92 



Wr[/, g] dxdx' 
1 



fix>,E)[gix<,E) 



fix^,E)g{0,E) 



x=x'=() 



Wr[/,g] 
fiO,E) 



f{0,E){g'{0,E) 



f{0,E)g{0,E) 
fiO,E) 
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In our case, f{x,E) — U{—E,x — S) where U is that parabohc cyhnder function [T] 
which decays for positive argument. Thus, 



PiE) 



U'i-E,-5) 
U{-E,-5) 



The Green's function is also of course the moment generating function. Thus, 
already at this stage, we can use the Green's function to recover the mean epidemic 
size. The answer is 



nP{n)dn 



Ar-3/2 



92 



dxdx' 



dne^''^^nK{x,x' ,n/N) 



x—x'—Q 



92 d 



N'/' ^P{E) 



x=x'=0;E=l/2 



E=l/2 



a=-l/2,x=-S 



We thus need an expression for U{a,x) for a near —1/2. We can get this by perturbing 
about the Gaussian solution at a — —1/2, writing 

where ^pl satisfies the inhomogeneous equation 

where 6a is the shift in a. Then, since the two modes of the homogenous equation are 
/i = and /2 = ^'^ dte* the solution for -0i that decays as cc — *■ +00 is: 

px />oo 



The derivative w.r.t. a is then 



da 



ln(;7(a,a:)) 



_ ^1 

a=-l/2 ^aV-O 



dx'e-'='"/^f2{x') - f2{x)e^'^^ / dx'e-""'"/ 



We can now differentiate w.r.t. x and get 

n = [-e-V4y2(^)_e-V2 / dx'e-^'"/^ + h{x)e--'"/^ 



,5-'/2 



1 + erf 
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reproducing of course our previous result, Eq. (|2])! 

The last order of business is to take the inverse Laplace transform. We have 



P{n) 



For the threshold case, d 



1 



0, we have 

U'{a,0) 
U{a,0) ' 



U{a,-6) 



The integral can then be done by residues, giving 

,-(2fe+l)n 



P{n) 



' OO 



1)! 



fc=0 



fc!2'= 



TT (1 



-2n 



)3/2 



in agreement with our previous result. For general S, however, one has to do the integral 
numerically. The most important information however, the asymptotic behavior for 
small and large n, can be gleaned analytically. 

We can get the small n expansion of the distribution function by using the large a 
expansion of the integrand. Using the Hankel Countour Integral: 



1 

27ri 



1 



T{z) 



we have 



P{n) 



i/2N _ 



1 



/2N 



7-fiC30 



dae° 



7 — 'iCXD 

-1 
rFi72) 

6. 



,-1/2 



-3/2 



3/2 



8r(i/2) 



-1/2 



1 

^16 
5 



128 



\n) 



Thus we have that the leading order behavior at small n is P{n) « (47™^)"^/^ 
independent of 5. One can also verify that this small n series reproduces the full 
P{n) for the (5 = case. 

The large-n asymptotics is clearly given by the ground state of the wall problem. In 
the limit 5 00, the wall, relative to the bottom of the potential, moves to —00 and 
the ground state energy goes to 1 /2 in our units. This translates to a decay rate of zero, 
once the exp{n/2N) factor is taken into account. As S decreases the wall moves closer 
to the potential minimum and the energy (and the decay rate) rise monotonically. The 
energy is 3/2 in our units when the wall hits the potential minimum at (5 = 0. This 
leads to the decay behavior exp{—n/N) at threshold, in accord with the full solution 
in this case. The energy continues to rise as (5 is decreased below 0, leading to a faster 
decay in n, with the energy diverging in the limit (5 — > — 00. 

For example, the second excited state of the harmonic oscillator has zeroes located 
a distance 1/\/2uj = 1 to the right and left of the energy minimum. Thus, for 5 = —1, 
we expect the decay exp{—2n/N). Since the normalized eigenvector in our units is 



=Af[x^- 1] e-"'/4 
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where J\f-^ = 2e"i/2 + y2^erfc(V2/2) so that {(f>'{l))^ = Uf^e-^/^ w 1.208, then for 
large n, 

P{n) - 1.208e-2«/^ 

For (5 — +1 on the other hand, the solution of the Schroedinger problem is given by the 
first zero of [/(a, —1) at ai = —0.88824 = leading to a decay exponent of 0.388824. 
The coefficient of the exponent is given by J7'(ai, — l)/(J^[/(ai, — 1)) = 0.4365. This 
can be seen in Fig. |10[ where the correct exponential fallofF in both cases is seen for 
large n. For n/N < 2, the effect of the higher eigenvectors leads to deviation from a 
pure exponential behavior. In general, the large n approximation is accurate as long 
as the next higher eigenvector has decayed. This "energy gap" is approximately 1 for 
large positive S, rises to 2 at 5 = and continues to rise as S becomes more negative. 
However as the ratio of the energy gap to the ground state energy fails as S decreases, 
the role of the higher excited states becomes more pronounced as S decreases. 

Examining the limit of large positive 6 in more detail, the energy is very slightly 
above 1/2, so there is an extremely small decay rate. The actual ground state energy 
can be calculated as follows: We first shift x by S, so that the quadratic potential is 
centered at the origin. Then, since the ground state energy in the presence of the wall 
is close to the wall-free value of 1/2, we can write 

where satisfies the inhomogeneous equation 

and e is the shift in the energy. As before, the two modes of the homogenous equation 
are /i = and /2 = er^ dte^ The solution for ipi that decays as x +oo 

is: 

^^^^eh{x) dx'Ae--'"/^h{x')-ef2{x) dx' Ae'^'" 



Jo Jx 

The second term dominates for large negative x, so that 

X 

and then the boundary condition that (5) = gives 



V27r 
Then, 

da de 6 

and 



^p'{-S) = A 



ASe-''/' 



2 V 2 

This gives for the leading order asymptotics of P(n) for large n, large S: 

S2 



V2Td^ \NV2tt J 
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so that the total probability of a major epidemic is the integral of P{n), which is 
Ri 1 — I/-R0, as expected, and 

in accord with our previous result. 

In the limit of large n and large negative 5, we get that the ground state energy 
is large, approximately E w 5^/4, so that a ~ —(5^/4. However, there are many 
states with approximately this energy. It is best to work directly from our integral 
formulation. For large negative 5, we can use a WKB ansatz to write U{a,x) ~ 
exp(\/^5(x/\/^) which yields 

^ = 5' = -y^T^ 

U{a,x) 

Expanding this in a power series for large a, we get 



S7 r(3/2) f6^\' 



Performing the integral over a gives 

^ ' fc!r(3/2 - A;)r(fc - 1/2) U / 

,n/2Ar -3/2 f - r(3/2) sin(7r(fc - 1/2)) / d^n 
~ ^ kin \4N 

k=0 ^ 
_ 6 -nS^/iN 



\/47rn3 

We see that we have successfully summed all the leading order contributions. For 
large negative i5, this is cut off at n's of order so we may drop the e"/^^ term. 
Then, we have the result for a pure constant drift, and so matches on to the subcritical 
distribution (which is the same as the SIR case). It of course reproduces the correct 
mean as well, since in the integral over n, small n's predominate, and 



f 

Jo 



'^''''^ 2V^n3/2 - i-S) 



This in turn matches on to the subcritical result, n « 1/(1 — Rq), as Rq approaches 
one from below. 

An interesting subtlety arise if we consider the zeroth moment of the distribution. 

In normal circumstances this would be unity, but the scaling behavior of P{n) dictates 
that the normalization integral is formally of order N~^^'^ and furthermore diverges, 
due to the n~3/2 behavior of P for small n. Nevertheless, if we blindly forge ahead, we 
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find 

Plot = / dnP{n) 

jQ 

= N-^/^P{l/2) 

,.-i/2 t/'(-l/2,-^) 
C/(-l/2,~<5) 

S 

Clearly this finite answer is the result of an analytic continuation. To understand its 
significance, let us consider the difference between P{n) and the distribution for the 
constant bias random walk with the same S. For small n these distributions as we have 
seen are identical, so the difference is integrable. Integrating the constant bias random 
walk, we get 

pCB ^ f°° 1 ^--nS^iN ^ 

*°* Jo Vlnn^' 2VN 

Thus, the difference is (5 + \S\)/{2y/N). This is zero for (5 < 0, which is correct, since 
even without the space-dependent drift toward the origin, every walker will eventually 
hit the origin. On the other hand, for 5 > the difference is S/^/N, which reflects the 
fact that with the added space-dependent drift all walkers are guaranteed to return to 
the origin, while without only a fraction I/Rq ~ 1 — do. Thus, looking at the 

difference between distributions provides an excellent way to make rigorous the concept 
of "major epidemics" , even slightly above threshold in the critical regime. Even below 
threshold, it highlights the added role of the space-dependent drift in reducing the 
probability of larger epidemics in favor of smaller ones. 

6. Concluding Discussion 

We have exhibited an exact expression for the mean epidemic size in the SIS model 
of endemic infection. We have evaluated this in the limit of large population size, 
and shown the crossover behavior that occurs in the vicinity of the critical infection 
number, Rq = 1. We have also calculated the distribution function for the epidemic 
size, again focussing on the crossover regime. 

It is important to note that the crossover behavior is universal, independent of the 
details of the model. What is important is the existence of two fixed points of the rate 
equation dynamics and a critical parameter where the two fixed points interchange 
stability. In the case the functional form of the mean epidemic size as a function of 
S, the scaled distance to the critical point will be the same, along with the scaling 
behavior with N. This is also true in general for any first passage time problem where 
the transition rates are constant (independent of N and location) at the transition. 
The first passage time, i.e. the mean (physical) time to extinction, in the SIS model 
also exhibits a crossover behavior at the transition, albeit different than that of the 
mean epidemic size. This is due to the fact that the transition rates in time are location 
dependent, p'^ = ak{N — k)/N, = f3k. As Doering, et al. |1] did not investigate the 
transition behavior of the mean extinction time, we for completeness present it here. 
The mean extinction time (starting with one infected) in the crossover regime is given 
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by 



oo ^ 

fJK 



2 Jo 

In particular, at threshold, the integral can be performed analytically at we get 

t{S = 0) « ^ (ln2A^ + 7) 

The SIR model, as has been demonstrated [THIHITU], exhibits a different scaling in 
the threshold regime since its fixed point structure is different. The SIR model has a 
line of fixed points at J = 0, but no fixed point at finite /. In fact, any tendency to 
immunity (or death, for that matter) will cause an otherwise SIS model to exhibit SIR 
behavior in the threshold regime for large enough N. This is due to the fact that the 
Time dependent bias in the SIR model, no matter how small in strength, overwhelms 
the space-dependent bias [10 for large enough N. 
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